LDA+Gutzwiller Method for Correlated Electron Systems: Formalism and Its 

Applications 
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We introduce in detail our newly developed ab initio LDA+Gutzwiller method, in which the 
Gutzwiller variational approach is naturally incorporated with the density functional theory (DFT) 
through the "Gutzwiller density functional theory (GDFT)" (which is a generalization of original 
Kohn-Sham formalism) . This method can be used for ground state determination of electron systems 
ranging from weakly correlated metal to strongly correlated insulators with long-range ordering. We 
will show that its quality for ground state is as high as that by dynamic mean field theory (DMFT), 
and yet it is computationally much cheaper. In additions, the method is fully variational, the 
charge-density self-consistency can be naturally achieved, and the quantities, such as total energy, 
linear response, can be accurately obtained similar to LDA-type calculations. Applications on 
several typical systems are presented, and the characteristic aspects of this new method are clarified. 
The obtained results using LDA+Gutzwiller are in better agreement with existing experiments, 
suggesting significant improvements over LDA or LDA+U. 

PACS numbers: 71.15.-m, 71.27+a, 71.15.Mb 



I. INTRODUCTION 

The density functional theory (DFT)[1, 2] is very suc- 
cessful in solid state physics and materials science. First- 
principles calculations based on this theory, using the lo- 
cal density approximation (LDA) or the generalized gra- 
dient approximation (GGA), have been well developed 
and widely accepted as a powerful theoretical tool for 
explaining and predicting ground state properties and 
electronic structures of a large amount of materials such 
as simple metals and band insulators. However both the 
LDA and GGA fail when they are applied to strongly cor- 
related electron systems, a very important class of mate- 
rials in condensed-matter physics. These materials con- 
tains unfilled d or f shells such as Cuprates, Manganites, 
Ruthenates, Fe-pnictides, Plutonium as well as the Heavy 
Fermion systems. In the last twenty years, many efforts 
have been made to improve the situation, new meth- 
ods, such as LDA+U[4], self-interaction corrected (SIC) 
LDA [5] and LDA plus dynamical mean field (DMFT) 
theory [6], have been proposed to provide new computa- 
tional tools for the quantitative study of the strongly cor- 
related materials. Those methods are quite successful in 
many aspects, nevertheless a method that is practically 
efficient and can capture the key feature of the correlation 
effect as well is still absent for the ground state studies. 

One of the main features in correlated electron sys- 
tems is that although the electrons in those narrow 3d or 
4f bands are delocalized they still show some atomic fea- 
tures, which manifest itself in the appearance of the Hub- 
bard band and the enhancement of the effective mass. In 
weakly correlated electron systems, electron states are 
delocalized in real space, exhibiting nearly free electron 
behavior which leads to good energy bands description. 
The derealization feature grants suitable electron den- 
sity dependent forms for correlation energy as presented 



in LDA and GGA since the electron distribution is not 
far from homogeneous electron gas. However, if electrons 
exhibit strong localization feature of the atomic orbitals, 
it is better to describe the electron states in real space. 
The presence of strong on-site correlations require proper 
treatment of atomic configurations, which is orbital de- 
pendent and plays important roles in determining the 
physical property in this case. Methods such as LDA+U 
[4] and LDA+DMFT [6] are proposed as remedies since 
this orbital-dependent feature is absent in both LDA and 
GGA. These methods start from similar Hamiltonians in- 
cluding on-site correlations but operate in different ways. 

In LDA+U method, the on-site interaction is treated 
in a static Hartree mean field manner, it is suited for 
strongly correlated systems with long-range ordering, 
such as the AF ordered insulators, but it fails for interme- 
diately correlated metallic systems. In DMFT method, 
the self energy which is purely local in space is obtained 
in a self-consistent way, which make the LDA+DMFT 
method the most accurate and reliable method now. 
However, the frequency dependent feature of the self en- 
ergy makes it very time consuming, and the full charge 
density self-consistency, which is very important for the 
accurate total energy calculation, is hard to be achieved. 

Looking back to the progress of analytical treatment 
of strongly correlated system, we can notice that the 
Gutzwiller variational approach (GVA) has been proofed 
to be quite efficient and accurate [7-9] for the ground 
state studies of many important phenomena, i.e. the 
Mott transition, ferromagnetism and superconductivity. 
This approach was first introduced by Gutzwiller to 
study the itinerant ferromagnetism in systems with par- 
tially filled d bands described by the Hubbard Model[10]. 
In this approach, a many body trial wave function was 
proposed, in which the weights of unfavorable atomic 
configurations are reduced according to the variational 
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parameters. Both itinerant and atomic features can be 
described spontaneously by this type of wave functions. 
Thus, an unified description from weakly to strongly cor- 
related system can be built up by the GVA, this grants 
its capability to accurately capture the essence of corre- 
lated systems. Various techniques have been developed 
to formulate this approach [7, 11-14] for different model 
Hamiltonians. The reliability and feasibility of GVA ap- 
plied to correlated systems have been demonstrated by 
these theoretical studies. 

In this article, we will show that the GVA can 
be naturally combined with the DFT. As the result, 
the LDA+Gutzwiller (simply called LDA+G hereafter) 
method [15] is proposed for practical calculations of cor- 
related electron systems. To understand the formal- 
ism, we will show that a generalized Gutzwillcr density 
functional theory (GDFT) can be established following 
the same spirit of Kohn-Sham (KS) formalism in the 
DFT. The GDFT itself is rigorous, however, its exchange- 
correlation functional is unknown. By introducing cer- 
tain approximation to the exchange-correlation energy in 
GDFT, the LDA+G method can be derived, very similar 
to the LDA or LDA+U methods derived by approxima- 
tion to the exchange-correlation term in the KS formal- 
ism. In order to show the validity and the advantage of 
this method, we will demonstrate that GVA is as accurate 
as DMFT for the ground state properties, but computa- 
tionally much cheaper. In addition, the present method 
is fully variational, which guarantees that many of the 
important physical quantities, such as the force or the 
linear response can be naturally obtained from the vari- 
ational principle. Detailed formalism of this method will 
be explicitly introduced here, and we will also show that 
a fully charge density self-consistent procedure can be 
carried out, which is quite crucial for the total energy cal- 
culations. Furthermore we will also show that LDA+G 
method is easy to be implemented into the existing codes, 
particularly if the LDA+U method is already available. 

This paper is organized as follows. In Sec. II, GVA 
is introduced for a multiband tight-binding Hamiltonian. 
Then we make detailed comparison between GVA and 
DMFT results in Sec. III. The combination of GVA with 
DFT and its derivation from GDFT will be presented 
in Sec. IV. In Sec. V, we apply the method to several 
typical systems and the results will be discussed. The 
proofs of some equations are put in the Appendix. 



II. GUTZWILLER VARIATIONAL APPROACH 

We start with the GVA for the ground state of corre- 
lated electron model systems. The detailed description of 
GVA has been presented by many authors, here we refer 
to reference [7] for the review. For generality, we consider 
a model system with a set of localized orbitals, such as d 
or / electrons, which can be described quite generally by 
the multiband Hubbard model. The Hamiltonian reads 
[14] 



H = H + H int = ]T t1>fciC ja ,+Y,Hi (1) 

and 

H i = ^ U°' a hiahi,,! (2) 
ct,<x'(cj+<x') 

where a is combined spin-orbit index of localized orbitals 
basis {4>o-} on site i, a = 1, . . . , 2TV (TV is orbital num- 
ber, e.g. TV — 5 for d electrons). The first part is just 
a tight-binding Hamiltonian extracted from LDA calcu- 
lation, and the second term is the local atomic on-site 
interaction in which only density-density correlations are 
taken into account for simplicity. For generalized on-site 
interactions, please refer to Ref [14]. 

We first examine the Hamiltonian in atomic limit (i.e, 
considering only the Hi term for single site). There are 
2TV different spin-orbitals and each spin-orbital could be 
either empty or occupied, thus totally 2 2N number of 
multi-orbital configurations Hi is diagonal in the 

space casted by all |T) configurations since the on-site 
interactions are density-density type. 

E iT = (T\Hi\T) = ]T Ur' (3) 
cr.cr'er 

is the interaction energy of configuration |T) for the i- 
th site. (For general interactions, the atomic part should 
be diagonalized, and the eigen vectors are linear combi- 
nation of |r)). Of course those possible configurations 
should not be equally weighted, and electrons tend to 
occupy those configurations which has relatively lower 
energy. For this purpose, we could construct projectors 
which project onto specified configurations |T) on site i 

m ir = \i,T)(i,T\ (4) 

with the normalization condition, 

^2 ™ir = 1 (5) 
r 

since all the configurations { |r) } form a locally complete 
set of basis. 

In Eq. (1), if the interactions are not presented, the 
ground state is exactly given by the Hartree uncorrelated 
wave function (HWF) l^o), which is a single determinant 
of single particle wave functions. However, after turning 
on the interaction terms, the HWF is no-longer an good 
approximation, since there are many energetically unfa- 
vorable configurations. In a physical view, to describe 
the ground state better, the weights of those unfavorable 
configurations should be suppressed. This is the main 
idea of Gutzwiller wave functions (GWF). GWF \^g) is 
constructed by acting a many-particle projection opera- 
tor on the uncorrelated HWF. 

|*G>=£|*0> 

73 =n^=nE A ^r (6) 

i i r 
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The role of projection operator P is to adjust the 
weight of each configuration through variational param- 
eters A ir (0 < A ir < 1)- The GWF falls back to non- 
interacting HWF if all Air = 1- On the other hand, 
if A^ = 0, the configuration T on site i will be totally 
removed. In this way, both the itinerant behavior of un- 
corrected wave functions and the localized behavior of 
atomic configurations can be described consistently, and 
the GWF will give a more reasonable physical picture of 
correlated systems than HWF does. 

The evaluation of GWF is a difficult task due to 
its many-body nature. There are many efforts in 
the literature, and the most famous one is Gutzwillcr 
approximation[10], introduced by Gutzwillcr along with 
his proposal of GWF. In this approximation, the inter- 
site correlation effect has been neglected and the physics 
meaning was discussed in Ref [7] and Ref [16]. The exact 
evaluation of the single-band GWF in one dimension[12] 
and in the limit of infinite dimensions [13] were carried 
out. It turns out that Gutzwiller approximation is exact 
in the latter case. Extensions to multi-band correlated 
systems using Gutzwiller approximation were carried out 
by J. Bunemann et al. [14]. Meanwhile Gutzwiller ap- 
proximation was proofed to be equivalent to slave-boson 
theories [17- 19] on a mean-field level for both one-band 
case [20] and multi-band case [21, 22]. 

The expectation value of Hamiltonian Eq. (1) in GWF 
is : 

(VE G |ff|* G ) = (g^o) m 
[ >G (* G |* G ) <*o|£ 2 |*o> 1 ' 

Using the Gutzwillcr approximation, in the limit of infi- 
nite dimensions, according to Ref [14] we have, 

+ E E iT m iT 

»,r 



(8) 



where is the weight of configuration T, 

m lV = (* G |m ir |*G> (9) 

and 

^/mr j mr'.Dp; r . 



E 



(10) 



with D£, r =< r'|C^|r >, < z la < 1. (See Appendix 
for details). 

In order to understand the above Gutzwillcr results 
properly, it is better to compare it with the Hartree-Fock 
scheme. For this purpose, here we give the Hartree-Fock 
expectation value of Hamiltonian (1) using HWF \^o), 

<#>o= E c'^^>')o+E( e - +Ae -)^+ c 

(ii) 



where C is a constant , and Ae^, which is proportional 
to interaction strength U , is a correction to the on-site 
energy (level shift) introduced by the static mean field 
treatment of the interaction term. 

Comparing Eq. (8) with (11), now it is clear that the 
main differences between the Gutzwiller and the Hartree 
approaches are: (1) There are orbital-related factors 
in the former associated with the hopping terms, which 
describe the renormalization of kinetic energy, while the 
kinetic energy in the Hartree approach is not renormal- 
ized; (2) The interaction energy in the Gutzwiller ap- 
proach is not simply scaled with the interaction strength 
U, but it is related to the configuration weights. While 
in the Hartree approach, the presence of interaction term 
will contribute simply to the on-site energy correction in 
proportional to U after the mean field treatment. 

The total energy under the GWF can be obtained by 
minimizing Eq. (8) with respect to configuration weights 
rriir, which now in fact are variational parameters. Since 
more variational parameters are presented in this ap- 
proach, the obtained ground state total energy is much 
better than that in HWF. In other words, by using the 
GWF, the obtained ground state total energy is further 
reduced due to the reduction of interaction energy, but 
in the cost of kinetic energy. The balance of two (gain 
and cost) is achieved by the energy minimization with 
respect to variational parameters. 

For the convenience of our following discussions, here 
we would like to generalize the formalism and make sev- 
eral definitions. Any operator A acting on the GWF, can 
be mapped to an corresponding Gutzwiller effective op- 
erator A G which acts on the HWF (rather than GWF), 
requiring that its expectation values is kept as the same, 

(* G |4* G ) - (p \V*AP\* ) = (*o|i G |*o) (12) 



here we have, 



P^AP 



(13) 



If the operator A is a single-particle operator, such as 



Ao = Aff ClC ia> (where A?f = {ia\A \ja')), 

then similar to the above procedure for the evaluation of 
kinetic energy in Eq. (8) , its Gutzwiller effective operator 
(in Gutzwiller approximation) can be written as, 



Aq — e a™ z ia zj a 'Cf a Cj a i + y^(i 



^CicrCiaAl" 

(14) 



where again Zi a are orbital dependent renormalization 
factors, which are determined through the configuration 
weights presented in the projector P. Here please note 
that the diagonal term and the hopping term should be 
treated separately (see the Appendix). 

Following the above general definition, we can now de- 
fine an Gutzwiller effective Hamiltonian H G which acts 
on HWF, 



HZ 



TTU 

n int 



(15) 
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such that the following equation holds, 

(* \H G \* ) - (*g\H\* g ) =E g (16) 

here the kinetic part H§ can be written out according to 
Eq. (14), and the interaction part is 



H 



■■it — / , ^iTfTliT 

i-r 



(17) 



for the density-density type interaction as discussed 
above. 

Now we are coming to a stage that we can solve the 
Gutzwiller problem easily through energy minimization. 
In practice, the minimization procedure will be done it- 
eratively with each loop being divided into two steps. 
The first step is to fix the Gutzwiller variational parame- 
ters mr and find the optimal HWF. As we know the H G 
for given mr, which is non- interacting, this step can be 
easily done by diagonalize it and fill the corresponding 
bands up to the Fermi level. Then in the next step, we 
will fix the HWF and optimize the energy respect to all 
the Gutzwiller variational parameters mr- The explicit 
equation can be written as: 



dE G 
dm 



i.r 



j,jjti crcr' 

dz 

dmi.v 



dmi 



(18) 



In this second step of calculations, for the lattice model 
with crystal periodicity, usually additional constrains can 
be adopted: (1) There is no site-dependency for z ia fac- 
tors and occupation number n i(7 , i.e, Zi a — Zj a and 
nia- — nj a ; (2) The charge on each orbital should be 
kept to be the same as that obtained by HWF (for pure 
density-density interaction as discussed in [14]), in other 
words we have: 

£ (r| CtC< CT |r) m,, r = n ta = n\ a = (n) (19) 



When all m^r are obtained, go back the first step to con- 
struct an new effective Gutzwiller Hamiltonian again. By 
this recursive method, all parameters rn^r and |\& ) can 
be obtained self-consistently. 

Typically, the second step of variations, i.e. the opti- 
mization of mir, is not so easy for multi-band systems, 
because a large number of non-linear equations need to 
be solved spontaneously. Fortunately, following the steps 
described in our previous publication [23], we are able 
to transfer the non-linear equations into linear equa- 
tions set, and furthermore a so called "adiabatic solution 
searching" procedure can be adopted. Those techniques 
will greatly reduce the computational cost and stabilize 
the calculations. 



III. COMPARISON OF GUTZWILLER 
APPROXIMATION WITH DYNAMICAL MEAN 
FIELD THEORY 

In this section, we will compare the results obtained by 
the GVA and that by the dynamical mean field theory 
(DMFT) for the single and two bands Hubbard model. 
With the careful comparison of kinetic energy, interac- 
tion energy and quasi-particle spectrum, we are going to 
clarify the following important issue: Can the Gutzwiller 
approximation (GA) capture the important "incoherent 
motion" of the correlated-electrons or not? This prob- 
lem is considered to be the biggest shortcoming of GA, 
which prevents it to be widely used in the first-principles 
calculations of strongly correlated materials. As we will 
show below, the GA can definitely capture the effect of 
"incoherent" motion in the ground state by its multi- 
configuration nature, which leads to very good agreement 
with the DMFT ground state results for both the kinetic 
and interaction energies. While the situation is not as 
good for the excited states, since the variational param- 
eters in the GA are determined by optimizing only the 
ground state energy, not that of excited states. There- 
fore, GA is a much better approximation for the ground 
state than for excited states. Within the frame of GA, it 
is difficult to construct the high energy excited states cor- 
responding to the upper and lower Hubbard bands. That 
is why in the green's function obtained by GA, we only 
have quasi-particle part and no Hubbard bands. While 
the problem only exists for the high energy excited states 
not for the ground state and the low energy quasi-particle 
states. 

We start from the multi-band Hubbard model (1), 
for the clarity we only keep the intra-orbital hoping 
tl'J = ti t j5 at(T i and neglect the on-site energy ej ;CT in the 
following (restoring them does not change the conclu- 
sions). To describe quasi-particle states, an important 
physical quantity is Z-factor. Actually, there are two 
different definitions of Z in literature . The first one 
is the renormalization factor of the effective band width 
for the quasi-particles, the second one is the weight of 
the coherent part in the electron green's function near 
Fermi surface. As we will show below, in GA the Z- 
factors obtained by the above two definitions match each 
other, while in DMFT they are quite different. In the 
following comparison, we compute Z within DMFT by 

Z = {1 — U^o)^ 1 ; which is quasi-particle weight. 

While in GA calculation we define Z — z 2 . 

Under Gutzwiller approximation, the quasi-particle 
states and quasi-hole states can be expressed as [24] 



\< H ) = 



* ) for Ska > MF 

VC ka \^ ) for Ska < Mf 



With the above trial wave function, the excitation energy 
can be calculated as 



5 



-O- Gutzwiller 
■ DMFT 
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FIG. 1: Comparisons of calculated Z-factor, total energy, dou- 
ble occupancy and kinetic energy of single band Hubbard 
model with half-filling. In GA the band renormalization fac- 
tor and quasi-particle weight coincident. Z-factor from DMFT 
denotes quasi-particle weight. Double occupancy is < n-|-nj > 



_i_ jpp/h _ 
±K ka - 



p/h 



#1$ 



P/ h \ 



(Ot h ) 



(20) 



for quasi-particle (quasi-hole) excitations. The above 
equation can be evaluated by GA as shown in Appendix, 
which leads to a simple expression of Green's function, 



lla 



IUJ - 



{Ska - M-f) 



ith 



FIG. 2: Comparisons of calculated Z-factor, total energy, dou- 
ble occupancy and kinetic energy of single band Hubbard 
model with occupation number n=0.9. 




lla = 



for e ka > fi F 

for £ ka < 



being the weight of the coherent part spectrum, which 
can also be evaluated to be equal to z 2 under GA as 
shown in Appendix. Therefore within GA the quasi- 
particle weight coincident with the renormalization fac- 
tor of the kinetic energy, thus dynamical informations are 
captured by variational approach. 

In Fig. 1-4, we compare the kinetic energy, interac- 
tion energy and Z -factor from GA and DMFT calcula- 
tions. We choose the non-interaction density of state to 
be p(e) — -^p\/D 2 — e 2 , which corresponds to Bethe lat- 
tice with infinite connectivity. Anderson impurity model 
in DMFT is solved by Lanczos method, which gives es- 
sentially exact results. As could be seen in Fig. 1^4, 
GVA captures ground state energies quite well for almost 
all correlation strength and band fillings. Although as a 
variational approach it targets at total energy and does 
not ensure the correctness of kinetic and interaction en- 
ergy in principle, we still observe quite good coincidence 



FIG. 3: Comparisons of calculated Z-factor, total energy, dou- 
ble occupancy and kinetic energy of single band Hubbard 
model with occupation number n=0.8. 



of kinetic and interaction energy respectively. We also 
notice that as the band degeneracy increase there is fur- 
ther coincident between GA and DMFT results as shown 
in the Fig. 4 for two band case. 

For half-filling n = 1.0 where there is Mott insulator 
transition for large U, the introduction of GA further ne- 
glects spatial correlation and under-estimate the absolute 
value of kinetic energy (Fig. 1). In Mott phase, GA gives 
vanishing double occupancy, while DMFT result always 
shows finite double occupancy due to spatial fluctuation. 
The flaw at large U could be traced back to the fact that 
the starting wave function of Gutzwiller projection is un- 
corrected Fermi liquid state l^o)- This fact is also part 
of the reasons for the important shortcoming of GA: it 
miscaptures the high energy excited states (seen from the 
overestimate of Z-factor in Fig. 1^4 for large U). For all 
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two band SU(N) interaction 




FIG. 4: Comparisons of calculated Z-factor, total energy, in- 
teraction energy and kinetic energy of two band Hubbard 
model with SU(N) interaction. 



interaction strengths and band fillings, Z-factor from GA 
is larger than DMFT treatment, i.e., the method gives 
more weight to low energy coherent part. Nevertheless, 
we recently demonstrated [25] that this overestimation of 
Z-factor can be further corrected by properly taking into 
account the contribution from excited states (namely the 
incoherent motion of electrons). 

To conclude this section, we make the assertion that 
as a cheap tool for correlated systems GA has fairly good 
energy resolution, particularly good for the ground state 
total energy with dynamical information included, but 
careful must be taken when dynamical information as- 
sociated with high energy excitations is trying to be ex- 
tracted from the GA results. 



IV. COMBINING DFT WITH GUTZWILLER 
VARIATIONAL APPROACH 

In this section, we will discuss how we can combine the 
DFT with the GVA. The discussions will be separated 
into three parts. The detailed formalism of LDA+G 
method are explicitly derived in the first part, which are 
used in realistic calculations. In the second part, we first 
introduce a general Gutzwillcr density functional theory 
(GDFT), and then we derive the LDA+G formalism from 
the firm base of GDFT. In such a way, we demonstrate 
the rigidity of this method. Finally, in the third part, 
we will discuss the on-site interactions and the double 
counting term. 



A. Formalism of LDA+Gutzwiller Method 

As we discussed above, the strong on-site correlation 
is underestimated in LDA. For those strongly correlated 



materials, in which the correlations play very important 
roles in determine the electronic structure, this underes- 
timation may lead to qualitative mistakes. One common 
procedure to overcome this problem is that we treat the 
interactions more explicitly on top of LDA level, just like 
what has been done in LDA+U or LDA+DMFT schemes. 
The starting effective Hamiltonian is usually written as: 



H 



H 



LDA 



H, 



dc 



(21) 



where Hlda is the LDA part Hamiltonian extracted from 
the standard LDA calculation, Hi nt is the on-site inter- 
action term, and Hd c is the double counting term rep- 
resenting the average orbital independent interaction en- 
ergy already included by LDA. 

Both LDA+U and LDA+DMFT methods start from 
the same Hamiltonian as shown above, however they 
treat the problem in different ways. In the LDA+U 
scheme, Hartree like mean field approximation is used 
to solve the above Hamiltonian, which can capture the 
orbital dependent physics (which is absent in LDA), but 
the dynamical correlation is still not included. While in 
the LDA+DMFT method, the purely local self energy 
is evaluated by solving an effective quantum impurity 
model mapped from the original lattice model. With 
the frequency dependent self energy, not only the ground 
state properties but also the dynamical response around 
the equilibrium can be considered by LDA+DMFT. Be- 
cause of the frequency dependency of the self energy, 
the LDA+DMFT is quite expensive in computational 
time. In many applications, we are only interested in 
the ground state properties and it is quite important to 
develop a new computational method for correlation ma- 
terials, which is as fast as LDA+U and can capture the 
dynamical correlation effect as well for the ground state. 

As we have proposed in the previous paper [15], an al- 
ternative way to solve the problem is to use Gutzwiller 
wave function rather than single determinant Hartree 
wave function. This approach is much cheaper than 
DMFT, but its quality is as good as DMFT for the 
ground state determination (as been shown in the last 
section), because it can capture the dynamical correla- 
tion effect due to the multi-configuration nature of the 
Gutzwiller wave function. More importantly, this ap- 
proach is fully variational, and can be easily combined 
with the DFT as will be discussed below. 

Now the goal is to solve the Hamiltonian (21) by the 
GVA. For this purpose, we need to discuss the Hamilto- 
nian in more detail. Since the problem to be addressed 
here is generally orbital-dependent, the effective Hamilto- 
nian should be written in a set of complete orbital basis, 
which are always available, such as wannier functions or 
atomic orbitals. These orbitals can be denoted by \ia), 
in which i is site index, a is spin-orbital index and C ia 
is the corresponding creation operator. 

Following the basic idea of LDA+U or LDA+DMFT 
approaches, the Hlda term in the effective Hamilto- 
nian (21) is regarded as single-particle operator, it can 
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be therefore expressed in terms of a complete set of or- 
bitals as 



HlDA — tij aa 'Cj a Cj a ' 
ij,aa' 

Ujaa' ={ia\H LDA \ja') 



(22) 



Suppose all the orbitals on the same site are correlated, 
and the interaction term can be written as 



Hm = E U 

iaa' (a^a') 



(23) 



in which fii a — Cj a Ci a . Now it is easy to see that the ef- 
fective Hamiltonian (21) has the same form as that shown 
in Hamiltonian (1). (Please note the double counting 
term only contribute to a constant uniform shift, and 
has no orbital-dependency, as will be addressed in the 
later part). Then following the steps discussed in the 
section I, we will be able to solve the problem. This is 
the scheme used in most of post-LDA techniques, where 
a tight-binding fit to LDA results are first obtained, and 
then local orbital dependent interaction terms are im- 
plemented and the problem with interaction should be 
solved by some many-body techniques. We can therefore 
call the above procedures as post-LDA plus Gutzwiller 
approaches, which has been recently used for several ex- 
amples [26-28]. 

However, our intention is to develop a complete 
LDA+G method with full charge self-consistency and 
without tight-binding fitting. Two important factors 
have to be considered for this purpose: (1) Realistic ma- 
terials consist of non-strongly-correlated bands as well, 
which can be treated nicely by LDA, and those strongly- 
correlated bands, which require the Gutzwiller step. 
Proper separation of two sets of energy bands is there- 
fore necessary; (2) Full charge density self-consistency 
need to be considered. These will be the main points for 
our following discussions. 

We first divide the complete orbital basis into localized 
and extended orbitals, and the interactions are added 
only for localized orbitals, for example d or / orbitals 
in transition metal or rare earth compounds . The local- 
ized and extended orbitals are labeled by { | ia) = Cj a \0)} 
and {\id) = Cj s \0)} respectively, and the completeness of 
orbital basis requires that: 



(24) 



Under the representation of this basis set, the Hlda 
reads 



H LDA = (J2\icr)(i<T\+J2\^)(^\) 



iS 



H LDA (J2 \ja')(ja'\+J2\jS')(j6'\ 

jo' jS' 



(25) 



As discussed in Section I, the GWF \ ^>g) is constructed 
from the HWF l^o) with proper projection. Any opera- 
tor acting on GWF can be mapped to an Gutzwiller effec- 
tive operator which acts on HWF instead. Since Hlda 
only consists of single particle operators, following the 
definition in Eq. (13), its corresponding Gutzwiller effec- 
tive Hamiltonian can be written as, 



HfuA = (E Zi l<r \ia) H + J2 \ lS )( lS \) 

ia iS 

JWE^^WI+EliW'l) (26) 

jo' jS' 

+ ^(1 - z1 a )\ia){ia\H LDA \i<j){ia\ 

io 

To derive this Gutzwiller effective Hamiltonian, it is es- 
sential to understand that for those non-interacting or- 
bitals, the corresponding renormalization fact Zis is equal 
to 1. This formula could be further simplified using the 
completeness condition (24) 

Hl DA = (E H + 1 - E HH) 

ia ia 

jo' jo' 

+ E^ 1 - z io)\ i(J ) (i(T\H LDA \ia)(ia\ 



and the interaction energy is given as, 

(* G |iJ mt |* G ) =^£, r m,r 



(27) 



(28) 



Now it is clear that the complete basis set defined at 
beginning is actually not necessary for realistic calcula- 
tions, because only the localized orbitals (i<r\ appear in 
the above equation. The interaction term are also de- 
fined only for those localized orbitals. We then come to 
a stage very similar to LDA+U, where localized orbitals 
are defined and interaction within those orbitals is sup- 
plemented. What is in additional to LDA+U scheme is 
that the kinetic energy of each local orbitals is renormal- 
ized by factor z ia which need to be determined in terms of 
configuration weights and configuration energy through 
the variational approach as shown below. 

In realistic calculations for solid crystals, it is more con- 
venient to carry out the calculations in reciprocal space, 
especially for those plane wave methods. The transfor- 
mation to the reciprocal space is quite straightforward, 
because the Gutzwiller approximation keeps the trans- 
lational symmetry unbroken. We first define the Bloch 
states of localized orbitals \ia) 



i^) = ^E e ^i^) 



(29) 
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Then the Gutzwiller effective Hamiltonian H*£ DA in k- 
space can be written as, 

Hf DA = (£ z a \ka)(ka\ + M(ka\) 

ka ka 

H LDA (J2 z a ,\k'a'){k'a'\ + l*V)<fcV|) 

fe'cr' k'a' 

+ Y,( 1 -4)\kv)(k'<r\H LDA \k'<T)(ka\ 

(30) 



kk'o 



Let's define the projector P = J2ka Pka = J2k a \k<r)(ka\ 
which projects onto the Bloch state of localized orbital, 
then the projection to the remaining delocalized orbitals 
is taken into account by 1 — P. For convenience, here 
we also define another projector Q = J2k <x z a\k<j){ka\. 
Then we have 

hSda = (1 - P + Q)Hlda{1 -P + Q) 

+ 2(1 - zl)\ka){k'a\H LDA \k'cj){kc7\ ( 31 ) 

kk' (7 

The total energy is obtained by evaluating the expec- 
tation value of the Hamiltonian, 

E{p) = (* |#£ DA |*o) + J2 Ermr - Edc 

r 

= (*o\{\-P + Q)H LDA {\-P + Q)\*o) (32) 



in which £1d A 



J2k( kcr \ H LD A \ka), n a = 
^2 k (^o\k(j}(ka l^o); and Ed c is the double count- 
ing energy. 

Now the remaining task is to minimize the total en- 
ergy functional with respect to variational parameters: 
uncorrelated wave function |*o) an d atomic configura- 
tion weight mr . Very similar to the familiar Kohn-Sham 
equation, the uncorrelated wave function |*o) in the pe- 
riodic lattice can be written as a simple Slater Determi- 
nant of single particle wave functions \ip n k}- Two sets of 
variational equations can be derived from minimization 
of Eq. (32) with respect to \tp n k) and mr, respectively. 
For the variation with respect to \ipnk), please note that 
the orbital occupation number n a — X^fc(V'nfc|-P/ccr|V ; nfc) 
also depends on \tp n k) inexplicitly. 



dE( P ) 



Z \ H LDA + 



dE{ P ) dz a 



dfnk(^nk\ dz a dn a 

=e n k\i>nk) 



Pka ~ HdcMnk) 



dE{ P ) 
dmr 



= E 



dE{p) dz a 
dz a d mr 



+ E r = 



(33) 



(34) 



in which 

~~Q~ =~(~2 X] fnk(lpnk\PkcrHLDA + HhDAPka^nk) 



nk 



(35) 



When deriving this equation, the following relation is 
used 



3(1 -P + Q)) 

--z a \ko){ko\ 

--(1-P + Q)\ka)(ka\ 

--{l-P + Q)P ka 



(36) 



There are several constraints. The wave functions should 
be orthogonal and normalized, the total configurations 
weight must be unity, and for pure density correlations 
the local densities will not be changed in GVA: 



(tpnk\i>n'k'} = $n,n'$k,k> 

r 

2 (r| C\C a |r> m r = {n„) G = (n„> 



(37) 



Through the above steps, we will be able to solve the 
problem for fixed LDA Hamiltonian Hld A - 

Now the question is how can we achieve self- 
consistency in the charge density. This step is very 
crucial, and the reason is the following. As we dis- 
cussed above, all electrons (both delocalized and local- 
ized) should be included in realistic calculations. How- 
ever those delocalized orbitals arc treated in LDA level, 
and localized states are treated by the LDA+G step. The 
modification of the localized state will in return affect the 
charge distribution of all other delocalized state, particu- 
larly the charge transfer process between the delocalized 
and the localized orbitals may happen. This is of course 
important physics. If it is not treated properly, differ- 
ent conclusions may be drawn, as already discussed in 
the example studies for Nai_ K Co02 [29], where several 
post-LDA plus DMFT studies give different results. 

The charge density self-consistency can be achieved 
easily as long as the charge density can be constructed, 
because the LDA Hamiltonian is determined by the elec- 
tron density. In the present LDA+G scheme, the elec- 
tron density can be constructed from the Gutzwiller wave 
functions by: 

p= (tt G |p|* G ) = (*o|p G |*o) (38) 
Since p is also a one-particle operator, similar to previous 
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steps, with the help of Eq. (14) we have 
P° = E z ita \ia)(ia\ + 1-^ l Z(T >H) 

i<7 icr 

k)(H(E ^AtfWA + 1 - E (39) 
+ E(! - z f<j)N (T )(* cr ll r )( r IN o ')(*' j l 



or we can write down the expression in the momentum 
space with the following simple expression, 



p G = (1-P + Q)\r)(r\(l-P + Q) 
+ J2(l-4)\kv)(k'<T\\r)(r\\k'a)(ka\ 



(40) 



kk 1 



Eq. (33), (34) together with Eq. (10) and (40) provide a 
self-consistent scheme, which is named LDA+Gutzwiller 
method by us. Eq. (33) is similar to the KS-equation 
in LDA or GGA, except that the Hamiltonian has been 
replaced by the corresponding effective Gutzwiller one 
with orbital-dependent terms. Eq. (34), which is used to 
determine the configuration weight, and Eq. (10), which 
determine the factors z a , are newly introduced by GVA. 

We illustrate our schematic self-consistent loops for 
LDA+G method in Fig. 5. Two main steps in this 
scheme are: (1) for fixed factor z a , solving the Kohn- 
Sham-like equation to get uncorrelated wave function 
^o- (2) for fixed \I/o> calculate the configuration weights 
rnr and then obtain factor z a . The iteration loops end 
when both electron density and renormalization factors 
are self-consistent. This scheme could be easily imple- 
mented in all kinds of existing ab initio codes, no matter 
what kinds of basis set are used for the wave- function, be- 
cause the essential computational requirement is just the 
calculations of projection to some local orbitals. If the 
LDA+U method is already available in the original code, 
the implementation of LDA+G method will be much eas- 
ier, since the local orbitals defined in LDA+U can be also 
used in the LDA+G formalism. 

After the charge density self consistency has been 
achieved, the ground state properties such as the sta- 
ble crystal structure, magnetic structure as well as the 
elastic properties can be calculated, which is quite simi- 
lar with the standard LDA procedure. Besides that, we 
can also obtain the density of states by LDA+G. Be- 
cause the band dispersion E n k obtained by LDA+G is 
for the quasi-particle excitations, we can derive two types 
of density of states (DOS) from the LDA+G band struc- 
ture. The first one is the quasi-particle DOS, which can 
be expressed as 



PQP M = E S ( U ~ E nk) 



(41) 



nk 



Experimentally the electronic part of the low tempera- 
ture specific heat is directly determined by the quasi- 
particle DOS and thus can be estimated by LDA+G 



Initial guscss of charge 
density and z-factors 
P{r),z a 



1_ 



construct Gutzwiller effec- 
tive Hamiltonian H^ DA 



solve Kohn-Sham-like equation 

\rrG i 9E(p) dz a p _ 
Hdc]\lpnk) = £nk\lpnk} 



T 



calculate electron density and 
local orbital occupancies p(r),n a 



T 



calculate dE/dz a 



obtain configuration 
weights by this equation 




calculate quantities: 
total energy ... 



FIG. 5: Flow chart of self-consistent loops for 
LDA+Gutzwiller method. 



for correlated materials. Another type of DOS is the 
integrated electronic spectral function, which is called 
electron DOS in LDA+G. Using z a obtained by the 
Gutzwiller approximation, the weight of the quasi- 
particle peak at E n k appearing in the electronic spectral 
function can be expressed as 



•-nk 



kn Q 



nk 



kr 



1 -P nk 



(42) 



Then the electron DOS can be obtained by the summa- 
tion over all k , which reads 



Pel M = E Z nk<5 (w - E nk ) 



(43) 



nk 



As we mentioned in the previous section, only the co- 
herent parts (quasi-particle peaks) can be captured by 
LDA+G not the incoherent parts (Hubbard bands). The 
electron DOS in LDA+G corresponds to the low energy 
part of the photo emission spectrum, which is mainly 
determined by the quasi-particle dynamics. 

In the end, we would like to comment on the rela- 
tionship between LDA+G method and LDA or LDA+U. 
First it can be easily figured out that, LDA+G method 
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falls back to DFT-LDA method spontaneously in non- 
correlated systems. This can be seen from the multi- 
band Hubbard model in which on-site interactions go to 
zero, then no configuration is energetically unfavorable 
and GWF falls back to HWF. In LDA+G method, this 
case means that the supplemented on-site interaction en- 
ergy and double-counting term are zero, then all local- 
ized orbitals are really delocalized and the corresponding 
renormalization parameters are Z a = 1. The formalism 
falls back exactly to DFT-LDA. 

LDA+Gutzwillcr method can also cover the LDA+U 
method when applying to strongly correlated insula- 
tors with long range ordering, in which LDA fails while 
LDA+U makes its success. Actually this is quite easy 
to be understood. In the strongly correlated insulator 
with long range ordering and integer occupation, such as 
the anti-ferromagnetic phase in the half filled Hubbard 
model, the unit cell is doubled by the AF order, which 
greatly reduces the local fluctuation among the atomic 
configurations and thus increases the z-factor to be close 
to unity. Again in the z—\ limit, the GWF returns back 
to HWF, the LDA+G energy functional is equivalent to 
that of LDA+U when the ^-factors approach one. 



B. Derivation from Gutzwiller Density Functional 
Theory 

In the last section, we have derived the LDA+G 
method in a physical but yet not rigorous way. In this 
section, however, we will derive the LDA+G method from 
a sound base. We will first introduce a Gutzwiller den- 
sity functional theory (GDFT), which is rigorous and ex- 
act, just like the Kohn-Sham (KS) formalism developed 
from the density functional theory (DFT). In the KS for- 
malism, as long as we know the functional of exchange- 
correlation energy E^fJ^ , we can solve the ground state 
problem. This is also true for GDFT, where we have E^ c 
in its rigorous form instead of E^ c s ' ■ Of course, exact E^ c 
is unknown, and certain kinds of approximation have to 
be used in realistic calculations. In the KS formalism, if 
the LDA is used for the E^ c s (« E^ A ), then the LDA- 
KS type formalism is realized; in addition, if LDA+U 
approximation is used for the E^ c s (~ E^ A+U ), the 
LDA+U method can be obtained. We will show here 
that the LDA+G method can be actually regarded as 
a " LDA + U approximation in GDFT formalism 1 ' , where 
LDA+U type approximation is used for the exchange cor- 
relation term. We will first establish the GDFT, then 
we derive the LDA+G formalism from the firm base of 
GDFT. 



1. DFT and Kohn-Sham 

It is helpful to recall the basis of DFT first. The 
Hohenberg-Kohn (HK) theorem(author?) [1] shows that 
the total energy of a interacting electron system can be 



defined as a universal functional in terms of electron den- 
sity p(r). The ground state energy is the global minimum 
of the functional. The electron density that minimizes 
the functional is the exact ground state electron density. 
The equations are written as 

E[p] = (*\H\*)=T[p} + E int [p] 

P = (*|p|*> (44) 

where is the ground state many-body wave function, 
T is the kinetic energy, and E int is the interaction en- 
ergy. At this stage for simplicity, the energy due to ex- 
ternal potential are not included in the formula (it will 
be supplemented later). This theorem is exact, but can 
not be used directly since the explicit form of this func- 
tional is unknown. This problem can be transferred to a 
equivalent Kohn-Sham (KS) problem by using the well- 
known Kohn-Sham ansatz (author?) [2], which is now 
become one of the most important basis of first princi- 
ple electronic structure calculations for solid states. In 
this ansatz, a reference system is introduced, whose ex- 
act Hamiltonian is still unknown, but we know that its 
ground state wave function can be exactly written as the 
HWF l^o)- (Therefore, the reference system here is ac- 
tually a non-interacting system because we know that its 
wave function is l^o))- As long as the charge density of 
the reference system p° matches the true ground state 
charge density (i.e, p° = p), then from the Hohenberg- 
Kohn DFT, the total energy of the true system can be 
reproduced through the reference system, 

E[p] =E KS lp°] = T KS lp°] + Eg V] + E™]A 

=<*o|T|* > + E« s {p } + ££V] (45) 
p=p° = (*o|p|*o> 

The KS kinetic energy T KS and the Hartrce energy 
E^ s of the reference system is different with the true 
kinetic and interaction energy, T and Ei nt , but impor- 
tantly their functional forms are known. KS's idea is 
simply to re-organize the total energy expression, such 
that those known parts can be treated explicitly and 
all unknown parts are moved into the third term called 
exchange-correlation energy E^f ■ It is in this sense that 
KS formalism is still exact for the ground state total en- 
ergy. The merit of KS formalism is that the problem is 
solvable as long as the functional form E^ c s is known. By 
definition, the exchange-correlation energy is given as, 

E™ = AT KS + AEfJ = (T — T KS ) + (E mt - < s ) 

(46) 

where two contributions should be physically included: 
(1) the correction to the kinetic energy AT KS ; (2) the 
correction to the interaction energy AE^f . 

The KS formalism up to now is exact, however certain 
kinds of approximations have to be made for E£ c in 
realistic calculations, as will be discussed below. Never- 
theless, here we want to put a note for the reference sys- 
tem. Once the approximation has been introduced for the 
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exchange-correlation term Eg , the nature of the refer- 
ence system may be modified. In reality, if the LDA+U 
approximation is used for the exchange-correlation po- 
tential in KS formalism, the reference system is no-longer 
non-interacting, and the l^o) is just a approximate wave 
function. 

2. Gutzwiller Density Functional Theory ( GDFT) 

Here we can establish a exact Gutzwiller density func- 
tional theory (GDFT) in parallel to KS formalism. One 
can see from KS formalism that, taking a non-interacting 
system as reference is not a necessity. What is really im- 
portant is to know the exact form of the wave-function of 
the reference system, such as \^o) in KS. The main ben- 
efit of this strategy is that in the reference system the 
kinetic operator can be evaluated explicitly in quite a 
simple form. In the spirit of Kohn-Sham ansatz, any sys- 
tems can be taken as reference providing that it has the 
same electron density as the true system. We can theoret- 
ically formulate the exact Gutzwiller density functional 
theory (GDFT), similar to KS, as follows. 

(1) To replace the original difficult interacting many- 
body system, we choose a auxiliary reference system, 
whose exact Hamiltonian is still unknown, but we know 
that its ground state wave function is given as the 
Gutzwiller wave function \^g) (rather than the HFW 
l^o))- An important point to be noticed here is that the 
fact whether the reference system is interacting or non- 
interacting actually does not matter, since the Gutzwiller 
wave function can be used to describe both the interact- 
ing and non-interacting systems in a better way. (The 
nature of the reference system depends on the choice 
of exchange-correlation potential as will be further dis- 
cussed in the following part); (2) Following the KS 
ansatz, we assume that the ground state density of the 
original interacting system is equal to that of the refer- 
ence system p G = p. (The representability is not rig- 
orously proofed at this stage, therefore this step is still 
an ansatz. However, considering the fact that the \t G 
automatically return back to \I/o in the non-interacting 
limit, it has the same spirit as KS ansatz which has been 
proofed to be valid for many applications.) (3) The ki- 
netic energy of the reference system can be explicitly 
written as T G = (* G |T|* G ). (4) All unknown parts are 
moved into the exchange-correlation energy Eg. Finally, 
the total energy and the charge density will be written 
as, 

E[ P ] =E G [p G ] = T G [p G ] + E G [p G ] + Eg[p G ] 

=<* G |f|* G ) + E G [p G ] + Eg[p G ] (47) 
p=p G = (^ G \p\H> G ) 

The above formulation of GDFT is still exact similar to 
KS. We can also come to the same conclusion as KS that 
if the exact exchange-correlation energy Eg is known, 
the exact ground state energy of the true system will 



be obtained. Again for the physical understanding, two 
terms are included in the exchange-correlation energy: 
(1) the correction to the kinetic energy AT G ; (2) the 
correction to the interaction energy AEf^ t , as expressed 
as, 

Eg = AT G + AE G t = (T — T G ) + (E mt - E G ) (48) 

As already mentioned, the *S?c automatically return 
back to in the non- interacting limit, therefore the 
present GDFT can be regarded as a general extension of 
original KS formalism. 

3. Approximations for E xc 

Certain approximation has to be introduced for the 
unknown exchange-correlation part in order to perform 
practically calculations. We will start from the LDA and 
LDA+U approximations used in KS formalism, and then 
we will show that by introducing the LDA+U type ap- 
proximation in GDFT, the LDA+G method can be de- 
rived. 

(1) LDA or GGA 

The most popularly used and widely accepted approx- 
imation is the LDA, where the exchange-correlation en- 
ergy is approximated as, 

Eg? « Eg DA = AT LDA + AEf n D t A (49) 

We assume the readers have the basic knowledge about 
LDA, we therefore do not discuss its details here. The 
only point we want to emphasize is that the LDA is ba- 
sically parametrized from the uniform electron gas, and 
only the local part of the exchange-correlation potential is 
kept, in other words, the non-local part of the exchange- 
correlation potential is neglected. This is the reason why 
LDA works well for simple metals, such as the Na, K, 
where wide s band cross the Fermi level, but it fails for 
strongly correlated systems, such as transition-metal ox- 
ides. 

(2) LDA + U for strongly correlated systems 

To overcome the problem of LDA for strongly cor- 
related systems, LDA+U method has been introduced, 
where the exchange-correlation energy is approximated 
as, 

E KS kE lda+u = AT LDA + AE LDA+U 

~AT LDA + AE^ A + (H int ) - E dc (50) 

— Exc A + (Hint) a — E dc 

The spirit of LDA+U method is that the LDA do 
not treat the interaction energy sufficiently well, and 
it need to be corrected for strongly correlated systems. 
Therefore, the interaction energy correction in original 
LDA AE^ A is replaced by the LDA+U counterpart 
^E^ A+U = AE\f t A + (JT jni ) - E dc . In such a way, 
the interaction term is treated more explicitly, and the 
energy is improved. 
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To understand the LDA+U formalism, it is very im- 
portant to notice the following points: (1) the reference 
system (in KS) is no longer non-interacting anymore, 
since the interaction term H int is supplemented from the 
exchange-correlation potential. Therefore, the \^o) is n0 
longer the rigorous eigen state of the reference system, 
instead it is just a approximation; (2) It is only up to 
this step that the definition of local orbital and interac- 
tion strength U is necessary. In such a way, the orbital- 
dependent potential is introduced (the LDA type E^ A 
only depends on the density, not on the orbital). 

The LDA+U method is quite successful for many of the 
insulating systems, which have AF long-range ordered 
ground state, however its quality is still not sufficient 
and can be improved further. The main drawbacks of 
LDA+U are two folds: (1) The supplemented interacting 
term Hi nt is treated by crude Hartree-scheme, and the 
interaction energy is typically overestimated. Actually, 
once the reference system becomes interacting system, 
the \&o is no-longer a rigorous wave function; (2) Due to 
the usage of \&o a s an approximate wave function, only 
the interaction energy is further corrected (over LDA), 
namely (Hi nt )o is the further correction to the interaction 
energy. But the kinetic part AT is still kept to be the 
same as that in LDA (AT LDA ). However it is known that 
the presence of interaction term should also renormalize 
the kinetic energy. Those drawbacks can be improved 
from our LDA+G formalism as will be derived below. 



4. Derivation of LDA+G from the GDFT 

The GDFT itself is also exact, however, certain ap- 
proximations have to be used for the E G C in practical 
calculations. Since the exchange-correlation energy is a 
functional of charge density, the easiest way of course is 
still to use the local density approximation, and neglect 
the non-local part of the potential, i.e, E G C « E^ A . For 
the Hartree energy, it only depends on the charge den- 
sity, and it is the same for both KS-DFT and GDFT, i.e. 
E% = E§ s . Therefore, after applying the LDA to E G C , 
all the potential energies in GDFT recover to be the same 
as that in LDA-KS. In this limit, we already know that 
the reference system is a non-interacting system, and the 
wave function \1/ G should return back to *o- Therefore, 
all the above GDFT formalism returns back to LDA-KS 
if the LDA is used for E G C . So far, we gain nothing from 
the usage of GDFT. However, if the LDA+U type ap- 
proximation is used for E G C in GDFT, the situation will 
be much improved. 

As have been discussed above, to overcome the prob- 
lem of LDA for strongly correlated system, the strategy 
of LDA+U approximation is to introduce a supplemented 
interaction term Hint m the exchange correlation poten- 
tial, such that the electron-electron interaction can be 
treated more explicitly beyond LDA. The reference sys- 
tem now is no longer non-interacting, but is still used 
to approximate the wave function of the reference sys- 



tem in the LDA+U KS formalism. The same approx- 
imation for the exchange correlation term can be also 
used in the GDFT, namely a interaction term Hi nt can 
be supplemented in the exchange correlation potential 
to describe the electrons in localized orbitals better (be- 
yond uniform electron gas in LDA). Again, the reference 
system is a interacting system now, however, what is dif- 
ferent in GDFT is that the wave function of reference 
system is given as ^ G rather than V&o- Of course, the 
GWF \& G is much better than HWF f° r interacting 
system, and it is in this sense that the formalism is im- 
proved over LDA+U. 

Therefore, using the similar LDA+U approximation in 
the GDFT, exchange-correlation energy and the corre- 
sponding total energy can be written as, 

E[p}=T G [p}+E H [p}+Eg[p} 

=(* G \f\Hf G )+E H \p]+E° c \p] 
Eg «£&? A+G - AT LDA + AE^ A+G (51) 

=AT LDA + AE LDA + {Hmt)a _ 
= Exc A + {Hint)G — Edc 

It differs from the LDA+U KS scheme in the follow- 
ing two points: (1) The supplemented interaction term 
is more precisely dealt with the Gutzwiller wave func- 
tion, namely the (H int )G is used instead of (H int )o; 
(2) Although the AT LDA is still used in the exchange- 
correlation functional, the kinetic energy is actually im- 
proved through the replacement of T KS by T G . The 
usage of AT LDA in E G C requires more discussions. It is 
known that the drawback of AT LDA is that it keeps only 
the local part and neglects the non-local part. There- 
fore, to make improvement, non-local correction should 
be supplemented. However, it is seen from our above 
GDFT formalism, the non-local correction to the kinetic 
energy has been naturally included through the replace- 
ment of (* |T|*o) by (tf G |T|* G ). Therefore, in the 
exchange-correlation part, the non-local correction to the 
kinetic energy is no-longer necessary, and only local-part 
need to be considered. This is the reason why AT LDA 
can be used for E G C . The usage of AT LDA in E G C can 
also guaranty that the present LDA+G formalism return 
back to the LDA+U solution in the static limit, where 
the z-factors approaching unity and \^f G ) approaching 
1*°}- 

Up to this stage, we have finished all the necessary 
steps, and we show that the original LDA+G formalism 
discussed in the last section can be derived from a more 
rigorous base, namely the "LDA+U type approximation 
for the exchange correlation potential in the exact GDFT 
formalism". Using the Eq. (51), the Hamiltonian (21) 
discussed in the last section under GWF will be recov- 
ered. To be more practical, here we will write down the 
final version of the necessary equations explicitly. The 
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(52) 



total energy and the exchange-correlations read, 

E[p] =<tt G |f|* G ) + E H [p] + Jv ext pd 3 r + Eg[p] 

Eg *Eg DA+G = Eg DA + (* G |ff jnt |* G > - E dc 
in which the electron density now is, 

p(r) = (* G |r><r|* G ) = (*o|p G |*o) (53) 
and the kinetic energy operator is 



where i is electron label. The Hartree interaction energy 
of electrons is 



Er -\S 



\r — r 



(55) 



If the supplemented interaction part of the Hamilto- 
nian is diagonal in configuration space (This is true if 
only density-density interactions are considered), the in- 
teraction term in the exchange-correlation energy reads 
according to Eq. (2), (3), (8), 

^> G \H int \^ G ) = (^g^H^g) =^^rm ir (56) 

i i;T 

For convenience, external potential energy, Hartree en- 
ergy and exchange-correlation energy (of LDA part) can 
be grouped together, and written as an functional form 
of charge density, 

E eHxc [ P ] = E H [p] + J V ext pd 3 r + Eg DA [p] (57) 

and the effective potential is defined as 

SE eHxc [p] 



V, 



eHxc 



bp 



(58) 



which is exactly the same as that in the LDA-KS formal- 
ism. 

Now the total energy is, 



E[p,m tT ] =<* G |T|tf G ) + E eHxc [p] +J2 E irmiT - E dc 

i-r 

= (*o|T G |* ) + EeHxM + Yj E a-mir - E dc 

(59) 

The total energy is a functional of uncorrelated wave 
functions |\&o) and the configuration weight mjr, and 
have just the same form as Eq. (32). Both should be 
variationally optimized as already shown in the last sec- 
tion. 

For a summary of this section, we have performed a 
explicit derivation just in the spirit of the Kohn-Sham 
ansatz to incorporating the Gutzwiller approach into 
density functional theory. With some reasonable approx- 
imations, we proof that this scheme gives the same result 
as what we developed in the last section within a simple 
physical interpretation. This provides a reliable founda- 
tion of this new LDA+G method. 



C. Interaction and Double Counting Terms 

Similar to all kinds of LDA+U or LDA+DMFT calcu- 
lations, the interaction and the double-counting terms 
still remain to be defined explicitly in the present 
LDA+G method, it is therefore also a kind of semi- 
empirical ab initio method in this sense. Nevertheless, 
we can in general follow the same definition used in 
LDA+U or LDA+DMFT method. From the physical 
point of view, we only consider the strong on-site inter- 
actions of localized orbitals. The interaction strength 
can be explicit expressed with Slater integrals (or called 
Slater-Condon parameters) (F°,F 2 ,---) in the atomic 
limit. However, in practice, it is more convenient to use 
"Kanamori parameters", U,U' ,J, and J', which are com- 
binations of Slater integrals. The general form of the 
on-site interactions can be written as [30] 

a a=£a' ,xx' 



J 



- V c f 



(60) 



iax Cia X C ia'x C i<*'x 



a=£a' ,x 



_ ^_ \ " t t 

where a denotes localized orbital and \ denotes spin. 
The first two terms are intra-orbital Coulomb interac- 
tion, and inter-orbital Coulomb interaction, respectively. 
The Hund's rule exchange coupling are divided into three 
parts: one is the longitudinal part (the third term) which 
only involves density-density coupling; the other two 
terms (the 4-th and 5-th terms) describe the spin flip 
and pair hopping processes respectively. In the atomic 
case, the relation U = U' + J + J' holds to retain the 
rotational invariancc in orbital space. And for typical 
d-orbital systems, where spin-orbital coupling is not so 
strong, the relation J — J' also holds, we therefore have 
U' = U — 2Jm general. 

In this article, we restricted ourselves to the pure 
density-density interactions for simplicity, and the inter- 
actions to be considered are 

U' 



Hi U ^ ^ TT'iaj l^ia J. + ^ ^ ] ^iax^ict' x' 
a a=£a' ,xx' 



J 



J2 n - 

a=£a' ,x 



(61) 



This on-site interaction Hamiltonian is already diagonal 
in the atomic configuration |T) space, and the corre- 
sponding configuration energy Er is a linear combination 
of U, U' and J [14]. 

With the on-site interactions determined, we come to 
the question how much of them arc taken into account 
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in LDA, that is, how to write down the double count- 
ing terms. It is known that these interactions goes into 
LDA in an average way without orbital dependence. As 
already discussed in LDA+U or LDA+DMFT methods, 
we just follow the common choice of the double counting 
terms as [4]: 

E d c[ni] = ^ —ni(rii - 1) 

1 J (62) 
-^77(n iT (n iT - l) + n il (n ll - 1)) 

i 

where rij is the total electron number of localized orbitals 
on the same site i, n, — + = ^2 ax ni ax , and U, J 
are spherically averaged interactions, which can be given 
as [6], 

U = / A U + 21U') (63) 



J=U-U' + J (64) 

where I is the angular momentum number of the corre- 
sponding localized shell. 

The Coulomb and exchange interactions, U and J, can 
in principle be obtained using corresponding Slater in- 
tegrals. However, in real materials, the bare electron- 
electron interactions must be screened, and the Slater 
integrals have to be renormalized. Therefore, it is a hard 
task to determine the effective U and J exactly. In prac- 
tice, usually two possible ways are followed: (1) deter- 
mining the parameters from available experimental infor- 
mation empirically; (2) calculating the parameters from 
constrained LDA method [4] and the linear response ap- 
proach [31]. Depending on different method, different 
values might be obtained, however the important strat- 
egy is that for single fixed parameter, the method should 
be able to explain all possible properties spontaneously 
and systematically, rather than using different parame- 
ters for different properties. It is only in this way, the 
obtained results can be justified. We should also no- 
tice that the interaction parameters also depend on the 
choice of local orbitals, because of the different screen- 
ing processes involved. For example, both atomic or- 
bitals and wannier functions can be used to define the 
local orbitals, however generally the effective interaction 
strength for atomic orbitals should be larger than that 
for wannier orbitals because the former is more localized. 
To construct the wannier orbitals, either the projected 
wannier method [32] or the maximally localized wannier 
function [33] can be used. 

Finally, a very similar approach has been recently pro- 
posed independently by K. M. Ho. et.al. to combine the 
Gutzwillcr approach with the DFT [34] . The spirit of our 
method and their proposal are almost the same, however, 
they differs in the definition of the interaction and double 
counting terms. It is still remained to be justified which 
way should be the best in future. 




U (eV) 



FIG. 6: The calculated Z-factor as function of U for SrV03 
using LDA+G method. 



V. RESULTS AND DISCUSSIONS FOR 
REALISTIC SYSTEMS 

The above proposed LDA+G method was imple- 
mented in our BSTATE (Beijing Simulational Tool for 
Atom Technology) code [35], which uses plane- wave 
ultra-soft pseudo-potential method, with projected wan- 
nier function for the definition of local orbitals. We have 
applied the method to study several typical systems, 
where the strong e-e interactions play important roles. 
They are non-magnetic metal SrV03, magnetic metal Fe 
and Ni, AF insulator NiO, and unconventional supercon- 
ductor Nai-zCoC^. Some of the results have been pub- 
lished [15, 29] with emphasis on particular issues in each 
example. On the other hand, the main purpose of this 
full paper is to present the whole formalism of LDA+G 
and to demonstrate its advantage, namely what knowl- 
edge can be gained beyond LDA or LDA+U. Therefore, 
to keep the completeness of our present paper, here we 
would like to concentrate on the physical consequence of 
our method by grouping all results together. We will dis- 
cuss all those results in a totally different manner, such 
that we can understand the LDA+G method better. 
1. Band-narrowing and mass renormalization 

An essential quantity included in the LDA+G method 
is the kinetic renormalization factor Z = z 2 a due to 
the dynamic correlation. The z a factors are orbital- 
dependent, and can be self-consistently obtained from 
the energy minimization. Under the Gutzwillcr approx- 
imation, the Z factor can be also understood as the 
quasi-particle weight. It has been widely recognized that 
LDA type calculations overestimate the band-width of 
correlated-electron systems. The error-bar could be as 
large as an order (such as in Heavy Fermion system) de- 
pending on the strength of correlation. We will show 
here, that this band-narrowing (or mass renormaliza- 
tion) physics can be correctly obtained from the LDA+G 
method. For example, SrV03 is a intermediately corre- 
lated metal with id-t\ g configuration. It has simple cu- 
bic perovskite crystal structure, and magnetic instabili- 
ties are not involved for the ground state property [36]. 



15 



00 

S 4 




Energy (eV) 



FIG. 7: The calculated density of states for SrVOa using LDA, 
LDA+U and LDA+G methods. We used the C/=5.0eV and 
J=1.0eV, and the electron DOS (rather than quasi-particle 
DOS) is shown, in the calculations of LDA+G. 



Although the LDA calculation can correctly predict the 
non-magnetic metallic nature of the ground state, the cal- 
culated band width is about 40% wider than photoemis- 
sion observation [37] , and the estimated effective mass is 
about 2-3 times lower than experimental results from spe- 
cific heat and susceptibility. On the other hand, all these 
features can be improved by LDA+G calculations, and 
correct band-narrowing and mass renormalization can be 
obtained as shown in Fig. 6 and Fig. 7 by using reason- 
able J7(~ 5.0eV). To gain further understanding, here we 
would also like to compare the results to that obtained 
by LDA+U method. In the LDA+U, only the interac- 
tion energy part is corrected over LDA, and the kinetic 
part is not renormalized, as already discussed in the pre- 
vious section. Therefore, LDA+U can not explain the 
observed large renormalization as shown in Fig. 7, where 
the density of states (DOS) obtained by LDA+U almost 
coincide with that by LDA. 

Finally, we will also show below that the band- 
narrowing and mass renormalization is such an impor- 
tant quantity and common advantage of LDA+G that it 
is encountered for all the examples we studied. 
2. Improved spin polarization 

Except the band renormalization, we will show here 
that the spin polarization of magnetic systems can be also 
improved to certain extend. To understand the physics 
better, we would like to divide the mechanisms of spin po- 
larization into two parts: (1) The inter-site exchange (or 
the spacial long-range exchange); (2) The intra-site ex- 
change (mostly the inter-orbital Hund's coupling). Such 
a separation is not rigorous, but just for the physical un- 
derstanding. It is important to note that, in our formal- 
ism, the spacial inter-site part still remain to be treated 
by LDA level, and only the intra-site part is improved 
explicitly. Of course, through the charge-density self- 
consistency, the inter-site part may be also tuned slightly, 
but it is not a main effect. It is therefore understood that 
the issues related to the inter-site exchange, such as the 
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FIG. 8: The calculated magnetic moment (per Fe) of bec FM 
Fe as function of U and J, by using different methods. 
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FIG. 9: The calculated DOS of bec FM Fe using different 
methods. The parameters U—7.0eV and J=1.0eV are used in 
these calculations. 



spin spacial fluctuation or the geometrically frustrated 
spin systems, can not be improved through the LDA+G 
treatment. 

Even for the intra-site interaction, it is treated both in 
LDA+U and in LDA+G, what will be the difference? In 
the LDA+U, it is treated from the static mean field level, 
which always tends to give larger spin polarization than 
that in LDA for positive effective U e ff(=U-J). On the 
other hand in the LDA+G method, the dynamic effects 
are included and the intra-site (inter-orbital) charge and 
spin fluctuations are all included in a better way. It is 
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in this sense that the results by LDA+G should be more 
reasonable. To demonstrate the effect of LDA+G on spin 
polarization, here we would like to show three examples: 
(1) For ferromagnetic (FM) bec Fe, as shown in Fig. 8 
and 9, the calculated magnetic moment in LDA+U is 
always larger than that in LDA, and is significantly 
overestimated compared to experiments. On the other 
hand, in LDA+G, this overestimation by LDA+U is sup- 
pressed. (2) For the bulk fee Ni, the calculated moment 
by LDA+G is even smaller than that in LDA, in bet- 
ter agreement with experiments (see Table I). (3) For 
Nai_a;Co02, the LDA level calculations predict that the 
system is magnetic for all the doping range, the LDA+U 
calculations even enhance the tendency to be magnetic, 
in contrast to experimental observation. However, using 
LDA+G, we show in Ref. [29] that the ground state is 
actually non-magnetic for the intermediate doping range 
(around 0.3< x <0.5), in nice agreement with experi- 
ments. 

3. Total energy and equilibrium properties 

A big advantage of LDA type calculations based on 
DFT is its ability to get the ground state total energy 
accurately. Here we will show that, by explicitly treat- 
ing the interaction term through our LDA+G formalism, 
the calculated total energy and equilibrium properties 
of correlated electron systems can be also improved sig- 
nificantly. Bulk Fc and Ni are typical magnetic met- 
als with intermediate correlations, where LDA produce 
big error bar for the ground state properties in compar- 
ison to experiments. For Fe, the LDA even fails to pre- 
dict the correct bee FM ground state (although GGA 
correctly do so, the reason is not clearly understood). 
The results, summarized in Table I shows that most of 
the discrepancies are systematically improved, compared 
with experiments, suggesting the advantages of present 
scheme. First of all, the bee FM ground state is now 
correctly predicted by LDA+G (see Ref. [15] for origi- 
nal figure), we therefore understand that the failure of 
LDA to predict the correct ground state is due to its un- 
derestimation of strong on-site correlation. Second, the 
calculated equilibrium volume, bulk modulus, magnetic 
moment, specific heat coefficient, and band width are all 
improved in a systematic way by a simple fixed interac- 
tion strength (U=7.0eV and J=1.0eV). This is in sharp 
contrast to that obtained in LDA+U, for instance, the 
LDA+U may also get the correct equilibrium volume by 
certain U value, but the obtained magnetic moment will 
be unreasonably larger than experimental results if the 
same U is used. 

4- Large gap AF ordered Insulator 

Now we come to discussions for the large gap AF or- 
dered insulator with integer occupation, where LDA+U 
works well. We will show that the LDA+G actually gives 
similar results in this limit. The reason is very straight 
forward as has been pointed out in the formalism. In the 
present LDA+G scheme, both the on-site level and the 
kinetic energy should be renormalized due to the presence 
of interaction term. However, in the case that long-range 



TABLE I: The calculated property parameters for bec FM 
Fe and fee FM Ni in comparison with experimental results. 
They are equilibrium lattice constant ao, bulk modulus B, 
spin magnetic moment M, specific heat coefficient 7, and the 
occupied energy band width W . The experimental data are 
from Ref. [38]. (This table is a reproduction of our results 
that published in Ref [15].) 







ao (bohrs) 


B(GPa) 


M( MB ) 




W(eV) 




LDA 


5.21 


227 


2.08 


2.25 


3.6 


Fe 


LDA+G 


5.39 


160 


2.30 


3.52 


3.2 




Exp. 


5.42 


168 


2.22 


3.1,3.69 


3.3 




LDA 


6.49 


250 


0.59 


4.53 


4.5 


Ni 


LDA+G 


6.61 


188 


0.50 


6.9 


3.2 




Exp. 


6.65 


186 


0.42,0.61 


7.02 


3.2 



ordering is established with integer occupation, if the en- 
ergy gap is big, each orbital should be close to cither fully 
occupied or totally empty, because the charge fluctuation 
between the states should be small. In this limit, the ki- 
netic renormalization Z factor will be very close to unity, 
and "Lq returns back to ^o- Therefore, the kinetic renor- 
malization is very small, and only the renormalization to 
the on-site level take effect, this is exactly just the limit 
that obtained in LDA+U. As we have shown in the cal- 
culations for NiO [15], the obtained electronic structure 
is very similar to that of LDA+U. However please note, 
even for the AF long-range ordered insulators, if the band 
gap is small and the spin moment is far away from in- 
teger, the dynamic processes crossing the band gap may 
also take effect, in this case, the Z factor will be no longer 
unity, and of course, the results by LDA+G will be differ- 
ent with that of LDA+U, and the one by LDA+G should 
be more close to reality. We expect that this situation 
may happen in the LaTiC>3 [39], where the gap is about 
0.2eV and the calculated moment from LDA+U is much 
larger than that observed experimentally. 
5. Effect of charge density self- consistency 

Here we will show that the charge density self- 
consistency is really important for the calculations on re- 
alistic systems, and we take Nai_ 2: Co02, a typical multi- 
orbital system, as an example. As we have pointed out 
in our recent publication [29], all the issues discussed 
above, such as the band-narrowing, spin polarizations, 
orbital fluctuations, are encountered in Nai_a;Co02, and 
systematic improvement are obtained through LDA+G 
treatment. However, we want to take Nai_ :c Co02 as an 
example to demonstrate the importance of charge self- 
consistency, because several post-LDA techniques (with- 
out charge self-consistency) have been applied to this 
compound and conflicting results are obtained [28, 40]. 
The issue is related to the relative splitting of energy level 
between e' g and a\ g states, and the appearance of e' g hole 
pockets at the Fermi surface (for x=0.3). If the splitting 
is large, the e' g orbital will be totally occupied, and there 
will be no e' g hole pockets at the Fermi surface. Starting 
from different Hamiltonians fitted to LDA band struc- 
ture, advanced techniques such as Gutzwiller or DMFT 
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FIG. 10: The calculated level splitting between the oi 9 and 
the e' g states for Nai-^CoO^. The original figure is obtained 
from Ref. [29]. The one- loop result corresponds to the calcu- 
lation without charge density self-consistency. 



have been applied in post-LDA scheme, however, one of 
the results suggests the absence of e' g hole pockets [28], 
and the other suggests the appearance [40] . It is now un- 
derstood [41] that the main reason is due to the difference 
in the fitted tight-binding Hamiltonian, namely the crys- 
tal field splitting (or the on-site energy) in the two studies 
are different. In our LDA+G method, full charge self- 
consistency is achieved and no tight-binding fitting is re- 
quired. Only after such kind treatment, the discrepancy 
can be nicely resolved [29]. On the other hand, because 
of the feedback effect in the charge self-consistency, the 
on-site level renormalization which is overestimated by 
post-LDA techniques is now suppressed as show Fig. 10. 

In summary, we have show in this full paper the de- 
tailed formalism of LDA+Gutzwiller method, and its 
firm derivation from the GDFT. By comparing the re- 
sults to that obtained by DMFT, we have shown that 
the energy resolution of Gutzwiller approach is pretty 
good for the ground state determination. It is computa- 
tionally cheaper, and yet with dynamic fluctuations in- 
cluded. The calculated results for several typical systems 
demonstrate that it can be widely applied to many of the 
correlated electron systems with quality beyond LDA+U. 
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APPENDIX 

In this appendix, we will proof several equations dis- 
cussed in the text part. We will pay special attention to 
when and how GA is applies. Even a variational ansatz 
is given for a lattice model, evaluation of expect value is 



not straightforward. GA is a systemically approximation 
to make the evaluation accessible. The spirit of GA ac- 
tually is to neglect Wick contractions between operators 
with different site/orbital indices, thus in the following 
we often need to Fourier transform the expression into 
real space and then apply GA. 



1. Proof of Eqn. 8 

First we note that by choose A^r = y^jiS I^g) 
is normalized under GA. (*g|*g) = n(*o|-F?l*o) = 
IIEr^<*oKr|*o> =n(£r™*r) = 1. In the first 

i zT i 

equality we separate the average of a projection operator 
string into the product of single site averages. 
Then the expectation value of kinetic energy is 

(* G |flo|*G> 



\ 



i j_ 



x (^o\m l . r 'C} a m ry r l m r x j C j(7 'm i . r ' |* > 

where in the first equality we adopt GA to neglect all 
Wick contractions from site i' ^ i,j and site i,j. To 
evaluate the expectation value we define 



"V; 

mi-Pi 
my 

J j j 



where l^T i ( lj-r •) are projection operators for orbital 
other than cr(cx')- Then we have 



J2 Jmr t m r ,m Tj m r , 



I, -A I, A 1,-A I, A 

to?. m°/m r m° , 
1 r. i r „ 



x (^o\niaCl (1 - n i(T ) (1 - n J<y .)C oal h oa ,\^ )D°, D°',* 
x (* \ClC ja ,\* ) 

= X;«i^*« '<*o|ctq,v|* ) (A.i) 
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The expectation value of interaction part is 



(H> G \H int \* G ) 



= EE^3 £ <*oK ;r |*o) 
i r 

= j2Y. E ? m ^ ( A -2) 



Put equations A.l and A. 2 together, we will have the 
equation (8) show in the text part. 



2. Proof of 7fc 



<7 **(J 



In this section, we proof jka = %a under GA, here we 
focus on quasi-particle sector without loss of generality. 
To compute < §\ a \ C\. a \ G >, we first Fourier transform 
the operators into real space then apply GA. 



= ^E e<fc(7 " J) ^ol C ^ C '/ t ^l*o) 



I, J 



I, J 

Ijtj 



+ |E eife(J " J) ^I^Ct CT |* ) 



i, j 

I=£J 



= z a (*o\C k „Cl\* ) 
= z a 6(s ka - Hf) 



Similarly, we could proof that under GA 

($jUC feCT |tf G ) = ^(*o|C| CT C fea |* ) = Z^iUF-Eka). 

Thus we have 7fc CT = z„ under GA. Also note that under 
GA the k dependence of the Z is missing, only the fact 
that above or below the Fermi surface matters. 



3. Evaluation of E v ha 



1 V 

N ^ 



z lkiI - J) E ^^<*o|^m J;r > /;rj Ctj* ) 



i,j 
ijtj 



T/ Tj 



- 1 ^(*o|C 7ct C] ct |* ) + ^ E ^ (/ - J) (*o|C Jff C] CT |* > 



A 7 



= 1 for £ feCT > ^L F 

Then the kinetic energy for spin a species reads, 



(^LlE^Q^L) 

= a7 E eife(W) ^t ii (*o\Cj„VCbC ja VCl\*o) 



i.j 



+ E eik[I ~ J) zZ^MCj.vctAvcl^^ 

J±I i,j 

For the first term since i ^= j, there is a constrain that 
I ^ i and I ^ j, other wise the expression vanishes. 
Then it equals to 



= (E 5 £ )4E^^oi^^^^i*o) 



For the second term, we have 4 cases: I = j but J ^ i; 
I 7^ j but J = i;I = j,J = i and / 7^ j, J ^ i. Following 
previous technique, one could find out each of the cases 
the projection operator V gives a z 2 a factor. Thus, 



< *LlE^ c ^l $ L> 

= zlY^^o^ClC^Cl^o) 
= 4/2 £ ^o\Cl a C kcT \^o)+zle ka 



First we show that the quasi-particle state is normal- 
ized under GA: 



While kinetic energy for \i 7^ a: 
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( $ LI Yl ^ C l C J^L) (*LI^/I*L> = Nj^Ermr + 0(1) 

r 

j i,j,n Put the kinetic and interaction energy together, we 

= [(*o|C / ^ / t JvI/ > + ^ e ^(/-./) ( v I , |C Ja C / t JvI/ >] get <<J>LI#I*L> = Ekn&Kl* < °\ c l c ^\° > 

j^j +N EY-rriY + z^eka + 0(1) , put the constant into chemical 

„ + potential, we have 

feM = z 2 a (e ka - Mf) 

Following the same routine, one could proof that For quasi-hole state, one could get similar results. 
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